{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbsphinx": "hidden"
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"import numpy as np\n",
"\n",
"from triqs.plot.mpl_interface import plt\n",
"plt.rcParams[\"figure.figsize\"] = (6,5) # set default size for all figures"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Square lattice susceptibility and the Random Phase Approximation (RPA)\n",
"\n",
"A nice example of two-particle response functions is the two-dimensional square lattice with nearest neighbour hopping $t$. Here we will use `TRIQS` and `TPRF` routines to compute the bare generalized susceptibility\n",
"\n",
"$$\\chi_{0, abcd}(i \\omega_n, \\mathbf{k}) ,$$ \n",
"\n",
"and the generalized susceptibility in the random phase approximation (RPA) \n",
"\n",
"$$\\chi^{(RPA)}_{abcd}(i\\omega_n, \\mathbf{k}) .$$\n",
"\n",
"In order to do this we first have to setup the tight binding Hamiltonian $H$, the single particle dispersion $\\epsilon(\\mathbf{k})$, and the single-particle Green's function $G_0(i\\omega_n, \\mathbf{k})$ of the system."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Tight binding Hamiltonian\n",
"\n",
"The tight binding Hamiltonian of the square lattice with nearest neigbhour hopping $t$ takes the form\n",
"\n",
"$$\n",
"H = -t \\sum_{\\langle i, j \\rangle, \\sigma} \n",
"\\Big( c^\\dagger_{\\sigma i} c_{\\sigma j} + \n",
" c^\\dagger_{\\sigma j} c_{\\sigma i} \\Big)\n",
"$$\n",
"\n",
"where $c^\\dagger_{\\sigma i}$ creates a fermion with spin $\\sigma \\in \\{\\uparrow, \\downarrow\\}$ on lattice site $i$, and $\\langle i, j \\rangle$ denotes summation over nearest neighbours $i$ and $j$.\n",
"\n",
"A representation of $H$ can be constructed using the `TBLattice` class in `triqs_tprf.tight_binding` where the hoppings are given as a dictionary with relative coordinate vectors as keys and hopping matrices as values. The unit vectors of the lattice, the position of the site local orbitals and names also needs to be setup, see below."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting on 1 Nodes at : 2019-06-05 17:03:01.636154\n"
]
}
],
"source": [
"from triqs_tprf.tight_binding import TBLattice\n",
"\n",
"t = 1.0\n",
"H = TBLattice(\n",
" units = [(1, 0, 0), (0, 1, 0)],\n",
" hopping = {\n",
" # nearest neighbour hopping -t\n",
" ( 0,+1): -t * np.eye(2),\n",
" ( 0,-1): -t * np.eye(2),\n",
" (+1, 0): -t * np.eye(2),\n",
" (-1, 0): -t * np.eye(2),\n",
" },\n",
" orbital_positions = [(0,0,0)]*2,\n",
" orbital_names = ['up', 'do'],\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Dispersion $\\epsilon(\\mathbf{k})$ and Fermi surface\n",
"\n",
"The Fourier transform of the real space hopping to momentum space $\\mathbf{k}$ gives the momentum space dispersion $\\epsilon(\\mathbf{k})$ and the alternative representation of $H$ as\n",
"\n",
"$$\n",
"H = \\sum_{\\mathbf{k}} \n",
"\\epsilon(\\mathbf{k}) c^\\dagger_{\\sigma \\mathbf{k}} c_{\\sigma \\mathbf{k}}.\n",
"$$\n",
"\n",
"To visualize $\\epsilon(\\mathbf{k})$, we can plot it along the high-symmetry path $\\Gamma - X - M - \\Gamma$ in the Brillouin zone. We use the helper function `k_space_path` to build the $\\mathbf{k}$-space path."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"kmesh = H.get_kmesh(n_k=(32, 32, 1))\n",
"\n",
"e_k = H.fourier(kmesh)\n",
"\n",
"G = [0.0, 0.0, 0.0]\n",
"X = [0.5, 0.0, 0.0]\n",
"M = [0.5, 0.5, 0.0]\n",
"\n",
"paths = [(G, X), (X, M), (M, G)]\n",
"\n",
"from triqs_tprf.lattice_utils import k_space_path\n",
"k_vecs, k_plot, k_ticks = k_space_path(paths, bz=H.bz)\n",
"\n",
"e_k_interp = np.vectorize(lambda k : e_k(k)[0, 0].real, signature='(n)->()')\n",
"\n",
"plt.plot(k_plot, e_k_interp(k_vecs))\n",
"plt.xticks(k_ticks, [r'$\\Gamma$',r'$X$',r'$M$',r'$\\Gamma$'])\n",
"plt.ylabel(r'$\\epsilon(\\mathbf{k})$')\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From this we can see that the bandwidth is $4t$ (with $t=1$) and that we have a Van Hove singularity at the $M$ point.To visualize the dispersion in the $(k_x, k_y)$ plane and the half-filled fermi surface where $\\epsilon(\\mathbf{k}) = 0$ we interpolate `e_k` on a two-dimensional grid and use `matplotlib.pyplot.contour` to plot the zero contour."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"e_k = H.fourier(H.get_kmesh(n_k=(32, 32, 1)))\n",
"k = np.linspace(-0.5, 0.5, num=200) * 2. * np.pi\n",
"Kx, Ky = np.meshgrid(k, k)\n",
"\n",
"e_k_interp = np.vectorize(lambda kx, ky : e_k([kx, ky, 0])[0, 0].real) \n",
"e_k_interp = e_k_interp(Kx, Ky)\n",
"\n",
"plt.figure()\n",
"extent = (k.min(), k.max(), k.min(), k.max())\n",
"plt.imshow(e_k_interp, cmap=plt.get_cmap('RdBu'),\n",
" extent=extent, origin='lower')\n",
"plt.colorbar()\n",
"plt.contour(Kx, Ky, e_k_interp, levels=[0])\n",
"plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I.e. the Fermi surface has a diamond shape with perfect nesting of the fermi surface for the momentum transfer vectors $\\mathbf{k}=\\pm(\\pi, \\pi)$ and $\\mathbf{k}=\\pm(\\pi, -\\pi)$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bare generalized susceptibility $\\chi_0$\n",
"\n",
"The bare generalized (particle-hole) susceptibility $\\chi_0$ is defined as the bubble diagram of a propagating particle hole excitation, i.e.\n",
"\n",
"$$\n",
"\\chi_0(\\tau, \\mathbf{r}) = - G(\\tau, \\mathbf{r}) G(-\\tau, -\\mathbf{r})\n",
"$$\n",
"\n",
"where $G$ is the single-particle Green's function.\n",
"\n",
"For a non-interacting system the single particle Green's function $G_0$ is known analytically\n",
"\n",
"$$\n",
"G_0(i\\omega_n, \\mathbf{k}) =\n",
"\\frac{1}{ i\\omega_n + \\mu - \\epsilon(\\mathbf{k}) }\n",
"$$\n",
"\n",
"and the product of two Green's functions can be performed analytically using the Lindhard formula\n",
"\n",
"$$\n",
" \\chi_0(\\omega_n, \\mathbf{q}) =\n",
" \\frac{1}{N_k} \\sum_{\\mathbf{k}} \n",
" \\frac{ f(\\epsilon(\\mathbf{k})) - f(\\epsilon(\\mathbf{k}+\\mathbf{q})) }\n",
"\t {i\\omega_n + \\epsilon(\\mathbf{k} + \\mathbf{q}) - \\epsilon(\\mathbf{k}) }\n",
"$$\n",
"\n",
"where $f(\\cdot)$ is the Fermi distribution function. \n",
"\n",
"For reference this formula is implemented in `TPRF`, for details on the generalization to multiorbital and accidentally degenerate dispersion see the documentation for `lindhard_chi00`."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from triqs_tprf.lattice import lindhard_chi00\n",
"from triqs.gf import MeshImFreq, Idx\n",
"\n",
"wmesh_lind = MeshImFreq(beta=5.0, S='Boson', n_max=1)\n",
"chi00_wk_analytic = lindhard_chi00(e_k=e_k, mesh=wmesh_lind, mu=0.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A computationally more efficient way of computing $\\chi_0$ is to perform the product of two single-particle Green's functions in imaginary time and real space and then fourier transform back to frequency and momentum space, i.e.\n",
"\n",
"$$\n",
"\\chi_0(i\\omega_n, \\mathbf{k}) = \n",
"\\mathcal{F}_{ \\{ \\tau, \\mathbf{r} \\} \\rightarrow \\{ i\\omega_n, \\mathbf{k} \\} } \n",
"\\{ -G(\\tau, \\mathbf{r}) G(-\\tau, -\\mathbf{r}) \\}\n",
"$$\n",
"\n",
"To do this we first compute the single particle Green's function of the lattice using the lattice Dyson equation\n",
"\n",
"$$\n",
"G(i\\omega_n, \\mathbf{k}) =\n",
"\\Big[ i\\omega_n + \\mu - \\epsilon(\\mathbf{k}) \\Big]^{-1}\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"from triqs_tprf.lattice import lattice_dyson_g0_wk\n",
"\n",
"wmesh = MeshImFreq(beta=5.0, S='Fermion', n_max=30)\n",
"g0_wk = lattice_dyson_g0_wk(mu=0., e_k=e_k, mesh=wmesh)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And then perform the bubble product using the dedictated function"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"╔╦╗╦═╗╦╔═╗ ╔═╗ ┌┬┐┌─┐┬─┐┌─┐\n",
" ║ ╠╦╝║║═╬╗╚═╗ │ ├─┘├┬┘├┤ \n",
" ╩ ╩╚═╩╚═╝╚╚═╝ ┴ ┴ ┴└─└ \n",
"Two-Particle Response Function tool-box \n",
"\n",
"beta = 5.0\n",
"nk = 1024\n",
"nw = 60\n",
"norb = 2\n",
"\n",
"Approx. Memory Utilization: 0.01 GB\n",
"\n",
"--> fourier_wk_to_wr\n",
"--> fourier_wr_to_tr\n",
"--> chi0_w0r_from_grt_PH (bubble in tau & r)\n",
"--> chi_wk_from_chi_wr (r->k)\n"
]
}
],
"source": [
"from triqs_tprf.lattice_utils import imtime_bubble_chi0_wk\n",
"chi00_wk = imtime_bubble_chi0_wk(g0_wk, nw=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Physical response function $\\chi_{S_z, S_z}$\n",
"\n",
"From the generalized susceptibility $\\chi_{abcd}(i\\omega_n, \\mathbf{k})$ we can obtain any two particle response function by contracting the rank 4 tensor indices with two quadratic operators. E.g. the spin-spin response along the z-axis is given by\n",
"\n",
"$$\n",
"\\chi_{S_z, S_z} = \\text{Tr}[ S_z \\, \\chi \\, S_z ]\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def chi_SzSz_contraction(chi):\n",
" \"\"\" Computes the trace Tr[ S_z \\chi S_z ]\"\"\"\n",
" Sz = np.diag([+0.5, -0.5])\n",
" chi_SzSz = chi[0, 0, 0, 0].copy()\n",
" chi_SzSz.data[:] = np.einsum('wqabcd,ab,cd->wq', chi.data, Sz, Sz)[:, :]\n",
" chi_SzSz = chi_SzSz[Idx(0), :]\n",
" return chi_SzSz"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To look at the structure of $\\chi_0$ we can plot it along the high symmetry path in the Brillouin zone (using the call operator of a triqs Green's function on `MeshBrZone`)."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def interpolate_chi(chi, k_vecs):\n",
" assert( k_vecs.shape[1] == 3 )\n",
" chi_interp = np.zeros(\n",
" [k_vecs.shape[0]] + list(chi.target_shape), dtype=complex)\n",
"\n",
" for kidx, (kx, ky, kz) in enumerate(k_vecs):\n",
" chi_interp[kidx] = chi((kx, ky, kz))\n",
"\n",
" return chi_interp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With some additional plot helper functions we can now look at the spin-spin response"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def plot_chi_1D(chi, label=None):\n",
"\n",
" chi_SzSz = chi_SzSz_contraction(chi)\n",
" chi_interp = interpolate_chi(chi_SzSz, 2*np.pi*k_vecs)\n",
"\n",
" plt.plot(ks, chi_interp.real, label=label)\n",
" \n",
" plt.grid()\n",
" plt.xticks(ticks=eps_k['K'], labels=[r'$\\Gamma$',r'$X$',r'$M$',r'$\\Gamma$'])\n",
" plt.title(r'Spin-response $\\chi_{S_z S_z}(\\mathbf{q}, \\omega=0)$')"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plot_chi_1D(chi00_wk, label=r'$U=0$')"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def plot_chi(chi, label=None):\n",
"\n",
" k = np.linspace(-0.75, 0.75, num=100) * 2.*np.pi\n",
" Kx, Ky = np.meshgrid(k, k)\n",
" k_vecs = np.vstack((Kx.flatten(), Ky.flatten(), 0*Kx.flatten())).T\n",
"\n",
" chi_SzSz = chi_SzSz_contraction(chi)\n",
" chi_interp = interpolate_chi(chi_SzSz, k_vecs)\n",
" chi_interp = chi_interp.real.reshape(Kx.shape)\n",
"\n",
" plt.imshow(\n",
" chi_interp, cmap=plt.get_cmap('magma'),\n",
" extent=(k.min(), k.max(), k.min(), k.max()),\n",
" origin='lower', vmin=0, vmax=chi_interp.max())\n",
"\n",
" plt.title(label); plt.colorbar()\n",
" plt.xlabel(r'$k_x$'); plt.ylabel(r'$k_y$')"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"label = r'Lindhardt response $\\chi^{00}_{S_z S_z}(\\mathbf{q}, \\omega=0)$'\n",
"plot_chi(chi00_wk, label=label)\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Random phase approximation (RPA)\n",
"\n",
"The full Bethe-Salpeter equation in the particle-hole channel has the form\n",
"\n",
"$$ \\chi^{(PH)} = \\chi_0 + \\chi_0 \\Gamma^{(PH)} \\chi^{(PH)} $$\n",
"\n",
"where all components are momentum and (three) frequency dependent quantities. In the random phase approximation (RPA) the vertex is apprximated by the (frequency independent) local interaction $\\Gamma^{(PH)} \\approx V$. \n",
"\n",
"The locality and frequency independence simplifies the BSE which formally is solved as\n",
"\n",
"$$ \\chi = [ 1 - \\chi_0 V ]^{-1} \\, \\chi_0 $$.\n",
"\n",
"In `TPRF`, there is the dedicated routine for RPA (`triqs_tprf.lattice.solve_rpa_PH`) that can treat the general multiorbital case.\n",
"\n",
"For generality, we take the nessecary steps to obtain the four index tensor $V$ from the operator expression of the Hubbard interaction."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"from triqs.operators import n, c\n",
"H_int = n(0, 0) * n(0, 1)\n",
"fundamental_operators = [c(0, 0), c(0, 1)]\n",
"\n",
"from triqs_tprf.OperatorUtils import quartic_tensor_from_operator\n",
"from triqs_tprf.OperatorUtils import quartic_permutation_symmetrize\n",
"V_int_abcd = quartic_tensor_from_operator(H_int, fundamental_operators)\n",
"V_int_abcd = quartic_permutation_symmetrize(V_int_abcd)\n",
"\n",
"V_abcd = np.zeros_like(V_int_abcd)\n",
"from itertools import product\n",
"for a, b, c, d in product(list(range(V_abcd.shape[0])), repeat=4):\n",
" V_abcd[a, b, c, d] = V_int_abcd[b, d, a, c]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can tune the interaction strength and compute the corresponding RPA generalized susceptibility."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"from triqs_tprf.lattice import solve_rpa_PH\n",
"\n",
"chi_wk_vec = []\n",
"U_vec = np.arange(1.0, 5.0, 1.0)\n",
"for u in U_vec:\n",
" chi_wk = solve_rpa_PH(chi00_wk, u * V_abcd)\n",
" chi_wk_vec.append(chi_wk)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plot_chi_1D(chi00_wk, label=r'$U=0$')\n",
"\n",
"for U, chi_wk in zip(U_vec, chi_wk_vec):\n",
" plot_chi_1D(chi_wk, label=r'$U=%2.2f$' % U)\n",
"\n",
"plt.legend();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the system has an instability at the $M$ point "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=(7, 3))\n",
"subp = [1, 2, 1]\n",
"\n",
"plt.subplot(*subp); subp[-1] += 1\n",
"label = r'Lindhardt response $\\chi^{00}_{S_z S_z}(\\mathbf{q}, \\omega=0)$'\n",
"plot_chi(chi00_wk, label=label)\n",
"\n",
"plt.subplot(*subp); subp[-1] += 1\n",
"label=r'RPA U=%1.0f $\\chi^{00}_{S_z S_z}(\\mathbf{q}, \\omega=0)$' % U_vec[-1]\n",
"plot_chi(chi_wk_vec[-1], label=label)\n",
"\n",
"plt.tight_layout()"
]
}
],
"metadata": {
"celltoolbar": "Edit Metadata",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.11"
}
},
"nbformat": 4,
"nbformat_minor": 2
}